homework 5, version 0

4.9 ms

Submission by: PoorRichRE (jazz@gmail.com)

xxxxxxxxxx
 
md"""
​
Submission by: **_$(student.name)_** ($(student.kerberos_id)@gmail.com)
"""

13.1 ms

6.1 ms

Homework 5: Epidemic modeling II

18.S191, fall 2020

This notebook contains built-in, live answer checks! In some exercises you will see a coloured box, which runs a test case on your code, and provides feedback based on the result. Simply edit the code, run it, and the check runs again.

For MIT students: there will also be some additional (secret) test cases that will be run as part of the grading process, and we will look at your notebook and write comments.

Feel free to ask questions!

1.9 ms

student

"PoorRichRE"

"jazz"

xxxxxxxxxx
 
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)
​
student = (name = "PoorRichRE", kerberos_id = "jazz")
​
# you might need to wait until all other cells in this notebook have completed running. 
# scroll around the page to see what's up

16.1 ms

Let's create a package environment:

23.6 μs

xxxxxxxxxx
 
begin
    using Pkg
    Pkg.activate(mktempdir())
end

263 ms

xxxxxxxxxx
 
begin
    Pkg.add(["PlutoUI", "Plots"])
​
    using Plots
    gr()
    using PlutoUI
    using Logging
    using StatsBase
end

359 s

In this problem set, we will look at a simple spatial agent-based epidemic model: agents can interact only with other agents that are nearby. (In the previous homework any agent could interact with any other, which is not realistic.)

A simple approach is to use discrete space: each agent lives in one cell of a square grid. For simplicity we will allow no more than one agent in each cell, but this requires some care to design the rules of the model to respect this.

We will adapt some functionality from the previous homework. You should copy and paste your code from that homework into this notebook.

58.4 μs

1.3 μs

Exercise 1: Wandering at random in 2D

In this exercise we will implement a random walk on a 2D lattice (grid). At each time step, a walker jumps to a neighbouring position at random (i.e. chosen with uniform probability from the available adjacent positions).

19.8 μs

Exercise 1.1

We define a struct type Coordinate that contains integers x and y.

21.9 μs

xxxxxxxxxx
 
struct Coordinate
    x::Int64 
    y::Int64
end

2.3 ms

👉 Construct a Coordinate located at the origin.

8.2 μs

origin

Coordinate x

xxxxxxxxxx
 
origin = Coordinate(0,0)

39.7 ms

Got it!

Good job!

6.5 ms

👉 Write a function make_tuple that takes an object of type Coordinate and returns the corresponding tuple (x, y). Boring, but useful later!

34.6 μs

make_tuple (generic function with 1 method)

xxxxxxxxxx
 
function make_tuple(c::Coordinate)
    return (c.x, c.y)
end

30.4 μs

Got it!

Great!

13.7 ms

Exercise 1.2

In Julia, operations like + and * are just functions, and they are treated like any other function in the language. The only special property you can use the infix notation: you can write

1 + 2

instead of

+(1, 2)

(There are lots of special 'infixable' function names that you can use for your own functions!)

When you call it with the prefix notation, it becomes clear that it really is 'just another function', with lots of predefined methods.

36.3 μs

xxxxxxxxxx
 
+(1, 2)

13.1 μs

+ (generic function with 206 methods)

xxxxxxxxxx
 
+

2.5 μs

Extending + in the wild

Because it is a function, we can add our own methods to it! This feature is super useful in general languages like Julia and Python, because it lets you use familiar syntax (a + b*c) on objects that are not necessarily numbers!

One example we've see before is the RGB type in Homework 1. You are able to do:
0.5 * RGB(0.1, 0.7, 0.6)
to multiply each color channel by $0.5$ . This is possible because Images.jl wrote a method:
*(::Real, ::AbstractRGB)::AbstractRGB

👉 Implement addition on two Coordinate structs by adding a method to Base.:+

21.6 μs

xxxxxxxxxx
 
function Base.:+(a::Coordinate, b::Coordinate)
    
    return Coordinate(a.x + b.x, a.y + b.y)
end

332 μs

Coordinate x

xxxxxxxxxx
 
Coordinate(3,4) + Coordinate(10,10) # uncomment to check + works

24.7 ms

Pluto has some trouble here, you need to manually re-run the cell above!

14.4 μs

Got it!

Good job!

64.0 ms

Exercise 1.3

In our model, agents will be able to walk in 4 directions: up, down, left and right. We can define these directions as Coordinates.

21.4 μs

possible_moves

Main.workspace3.Coordinate1Coordinate x

2Coordinate x

3Coordinate x

-1

4Coordinate x

-1

xxxxxxxxxx
 
# uncomment this:
​
possible_moves = [
    Coordinate( 1, 0), 
    Coordinate( 0, 1), 
    Coordinate(-1, 0), 
    Coordinate( 0,-1),
 ]

64.3 ms

👉 rand(possible_moves) gives a random possible move. Add this to the coordinate Coordinate(4,5) and see that it moves to a valid neighbor.

9.1 μs

Coordinate x

xxxxxxxxxx
 
Coordinate(4,5)+rand(possible_moves)

62.9 ms

We are able to make a Coordinate perform one random step, by adding a move to it. Great!

👉 Write a function trajectory that calculates a trajectory of a Wanderer w when performing n steps., i.e. the sequence of positions that the walker finds itself in.

Possible steps:

Use rand(possible_moves, n) to generate a vector of n random moves. Each possible move will be equally likely.
To compute the trajectory you can use either of the following two approaches:
1. 🆒 Use the function accumulate (see the live docs for accumulate). Use + as the function passed to accumulate and the w as the starting value (init keyword argument).
2. Use a for loop calling +.

7.3 ms

Main.workspace3.trajectory

xxxxxxxxxx
 
"""Retuns resulting Coordinate after n::Int random steps, starting from w::Coordinate initial position"""
function trajectory(w::Coordinate, n::Int)
    return accumulate(+, rand(possible_moves, n); init=w)
end

25.2 ms

test_trajectory

Main.workspace3.Coordinate1Coordinate x

2Coordinate x

xxxxxxxxxx
 
test_trajectory = trajectory(Coordinate(4,4), 2) # uncomment to test

985 ms

11.4 s

Got it!

Awesome!

278 ms

plot_trajectory! (generic function with 1 method)

xxxxxxxxxx
 
function plot_trajectory!(p::Plots.Plot, trajectory::Vector; kwargs...)
    plot!(p, make_tuple.(trajectory); 
        label=nothing, 
        linewidth=2, 
        linealpha=LinRange(1.0, 0.2, length(trajectory)),
        kwargs...)
end

197 μs

xxxxxxxxxx
 
let
        long_trajectory = trajectory(Coordinate(4,4), 1000)
​
        p = plot(ratio=1)
        plot_trajectory!(p, long_trajectory)
        p
    end
​
# ^ uncomment to visualize a trajectory

535 ms

Exercise 1.4

👉 Plot 10 trajectories of length 1000 on a single figure, all starting at the origin. Use the function plot_trajectory! as demonstrated above.

Remember from last week that you can compose plots like this:

let
    # Create a new plot with aspect ratio 1:1
    p = plot(ratio=1)

    plot_trajectory!(p, test_trajectory)      # plot one trajectory
    plot_trajectory!(p, another_trajectory)   # plot the second one
    ...

    p
end

18.3 μs

xxxxxxxxxx
 
let
    p = plot(ratio=1)
    steps = 1000
    starting = Coordinate(0,0)
    for i in 1:10
        random_walk = trajectory(starting, steps, 30)
        plot_trajectory!(p, random_walk)
    end
    p
end

172 ms

Main.workspace3.Coordinate1Coordinate x

2Coordinate x

3Coordinate x

-1

4Coordinate x

-1

-1

5Coordinate x

-1

6Coordinate x

7Coordinate x

8Coordinate x

9Coordinate x

10Coordinate x

xxxxxxxxxx
 
trajectory(Coordinate(0,0), 10, 30)

90.5 μs

Exercise 1.5

Agents live in a box of side length $2 L$ , centered at the origin. We need to decide (i.e. model) what happens when they reach the walls of the box (boundaries), in other words what kind of boundary conditions to use.

One relatively simple boundary condition is a collision boundary:

Each wall of the box is a wall, modelled using "collision": if the walker tries to jump beyond the wall, it ends up at the position inside the box that is closest to the goal.

👉 Write a function collide_boundary which takes a Coordinate c and a size $L$ , and returns a new coordinate that lies inside the box (i.e. $[- L, L] \times [- L, L]$ ), but is closest to c. This is similar to extend_mat from Homework 1.

27.9 μs

collide_boundary (generic function with 1 method)

xxxxxxxxxx
 
 function collide_boundary(c::Coordinate, L::Number)
    x0, y0 = 0,0  # placeholders for payload
    x0 = abs(c.x) <= L ? c.x : L  # if c is inside the LxL box
    x0 = (c.x < 0) & (x0 == L) ? -x0 : x0  # if c is outside the LxL box, handle signs
    y0 = abs(c.y) <= L ? c.y : L
    y0 = (c.y < 0) & (y0 == L) ? -y0 : y0
    return Coordinate(x0, y0)
 end

90.8 μs

Coordinate x

xxxxxxxxxx
 
collide_boundary(Coordinate(3,2), 10) # uncomment to test

46.9 ms

Exercise 1.6

👉 Implement a 3-argument method of trajectory where the third argument is a size. The trajectory returned should be within the boundary (use reflect_boundary from above). You can still use accumulate with an anonymous function that makes a move and then reflects the resulting coordinate, or use a for loop.

38.5 μs

Main.workspace3.trajectory

xxxxxxxxxx
 
​
"Retuns resulting Coordinate after n::Int random steps, starting from w::Coordinate initial position within LxL boundaries "
function trajectory(c::Coordinate, n::Int, L::Number)
    return [collide_boundary(x, L) for x in accumulate(+, rand(possible_moves, n); init=c)]
 end

148 μs

1.4 μs

Exercise 2: Wanderering Agents

In this exercise we will create Agents which have a location as well as some infection state information.

Let's define a type Agent. Agent contains a position (of type Coordinate), as well as a status of type InfectionStatus (as in Homework 4).)

(For simplicity we will not use a num_infected field, but feel free to do so!)

13.1 μs

xxxxxxxxxx
 
@enum InfectionStatus S I R

36.6 ms

x
 
abstract type AbstractAgent end

7.4 μs

Agent

x
 
# define agent struct here:
begin
    mutable struct Agent <: AbstractAgent
        status::InfectionStatus
        num_infected::Int64
        position::Coordinate
        path::Array  # coordinate history
    end
​
    Agent() = Agent(S, 0, Coordinate(0,0), [Coordinate(0,0)]) # default starting position is origin
    Agent(s::InfectionStatus, n::Int64, p::Coordinate) = Agent(s,n,p,[p]) # starting point will always be the first coordinate in path history
end

2.6 ms

Exercise 2.1

👉 Write a function initialize that takes parameters $N$ and $L$ , where $N$ is the number of agents abd $2 L$ is the side length of the square box where the agents live.

It returns a Vector of N randomly generated Agents. Their coordinates are randomly sampled in the $[- L, L] \times [- L, L]$ box, and the agents are all susceptible, except one, chosen at random, which is infectious.

33.7 μs

initialize (generic function with 1 method)

x
 
 function initialize(N::Number, L::Number)
    agents = [Agent(S, 0, Coordinate(rand(-L:L),rand(-L:L))) for i in 1:N]
    rand(agents).status = I
    return agents
 end

100 μs

Main.workspace31.Agent1Agent status

I::InfectionStatus = 1

num_infected

positionCoordinate x

pathMain.workspace3.Coordinate1Coordinate x

2Agent status

S::InfectionStatus = 0

num_infected

positionCoordinate x

-5

-5

pathMain.workspace3.Coordinate1Coordinate x

-5

-5

3Agent status

S::InfectionStatus = 0

num_infected

positionCoordinate x

-8

-1

pathMain.workspace3.Coordinate1Coordinate x

-8

-1

xxxxxxxxxx
 
initialize(3, 10)

99.7 ms

Got it!

Great!

349 ms

color (generic function with 1 method)

xxxxxxxxxx
 
# Color based on infection status
color(s::InfectionStatus) = if s == S
    "blue"
elseif s == I
    "red"
else
    "green"
end

66.9 μs

position (generic function with 2 methods)

x
position(a::AbstractAgent) = last(a.path) #a.position # uncomment this line

38.0 μs

color (generic function with 3 methods)

xxxxxxxxxx
 
color(a::AbstractAgent) = color(a.status) # uncomment this line

445 μs

Exercise 2.2

👉 Write a function visualize that takes in a collection of agents as argument, and the box size L. It should plot a point for each agent at its location, coloured according to its status.

You can use the keyword argument c=color.(agents) inside your call to the plotting function make the point colors correspond to the infection statuses. Don't forget to use ratio=1.

9.2 μs

visualize (generic function with 1 method)

xxxxxxxxxx
 
 function visualize(agents::Vector, L)
    scatter([c.x for c in position.(agents)],[c.y for c in position.(agents)], c=color.(agents), ratio=1)
 end

147 μs

xxxxxxxxxx
 
let
    N = 20
    L = 10
    visualize(initialize(N, L), L) # uncomment this line!
end

789 ms

Exercise 3: Spatial epidemic model – Dynamics

Last week we wrote a function interact! that takes two agents, agent and source, and an infection of type InfectionRecovery, which models the interaction between two agent, and possibly modifies agent with a new status.

This week, we define a new infection type, CollisionInfectionRecovery, and a new method that is the same as last week, except it only infects agent if agents.position==source.position.

10.0 μs

x
 
abstract type AbstractInfection end

12.9 μs

xxxxxxxxxx
 
struct CollisionInfectionRecovery <: AbstractInfection
    p_infection::Float64
    p_recovery::Float64
end

2.2 ms

Write a function interact! that takes two Agents and a CollisionInfectionRecovery, and:

If the agents are at the same spot, causes a susceptible agent to communicate the desease from an infectious one with the correct probability.
if the first agent is infectious, it recovers with some probability

20.8 μs

interact! (generic function with 1 method)

x
 
function interact!(agent::Agent, source::Agent, infection::CollisionInfectionRecovery)
    if agent.position == source.position
        # agents must be in the same location to interact
        if agent.status==S && source.status==I && rand() <= infection.p_infection
            # for potential infect, one must be S and the other I 
            agent.status = I
            source.num_infected += 1
        end
    elseif (agent.status == I)  && (rand() <= infection.p_recovery)
            agent.status = R 
    end
end

78.7 μs

1Main.workspace31.Agent1Agent status

S::InfectionStatus = 0

num_infected

positionCoordinate x

pathMain.workspace3.Coordinate1Coordinate x

2Agent status

I::InfectionStatus = 1

num_infected

positionCoordinate x

pathMain.workspace3.Coordinate1Coordinate x

2Coordinate x

xxxxxxxxxx
 
let  # test interact! function 
    agents = [Agent(S, 0, Coordinate(0,1)), Agent(I, 0, Coordinate(0,1))]
    interact!(agents[1], agents[2], CollisionInfectionRecovery(.5, .5))
    agents, position(agents[1])
end

244 μs

Exercise 3.1

Your turn!

👉 Write a function step! that takes a vector of Agents, a box size L and an infection. This that does one step of the dynamics on a vector of agents.

Choose an Agent source at random.
Move the source one step, and use collide_boundary to ensure that our agent stays within the box.
For all other agents, call interact!(other_agent, source, infection).
return the array agents again.

19.4 μs

1Agent status

S::InfectionStatus = 0

num_infected

positionCoordinate x

pathMain.workspace3.Coordinate1Coordinate x

2Coordinate x

-1

3Coordinate x

4Coordinate x

5Coordinate x

6Coordinate x

7Coordinate x

8Coordinate x

9Coordinate x

10Coordinate x

11Coordinate x

2Coordinate x

xxxxxxxxxx
 
begin  # get the ending position after 10 steps in L=5 
    steps = 10
    L = 5
    agent=Agent()
    agent.path = append!(agent.path, trajectory(agent.position, steps, L))
    agent, position(agent)
end

1.1 ms

rand_move! (generic function with 1 method)

x
 
function rand_move!(agent::AbstractAgent, steps::Number, L::Number)
"""move agent in a random direction of steps::Number, within L::Number boundaries"""
    agent.path = append!(agent.path, trajectory(agent.position, steps, L))
    agent.position = position(agent)
end

60.6 μs

Agent status

S::InfectionStatus = 0

num_infected

positionCoordinate x

-1

pathMain.workspace3.Coordinate1Coordinate x

2Coordinate x

-1

x
 
let  # test rand_move! function
a = Agent()
rand_move!(a, 1, 5)
a
end

45.4 ms

step! (generic function with 1 method)

xxxxxxxxxx
 
 function step!(agents::Vector, L::Number, infection::AbstractInfection)
    # choose a source at random
    n_agents = length(agents)
    source_idx = rand(1:n_agents)
    source = agents[source_idx]
    # move the source by step
    step = 1
    #@debug "source=agent $(source_idx), start-$(source.position)"
    rand_move!(source, step, L)
    #@debug "end-$(source.position)"
    # interact source with all other agents
    for i in 1:n_agents
        if i != source_idx
            interact!(agents[i], source, infection)
        end
    end
    return agents
 end

80.3 μs

Main.workspace31.Agent1Agent status

S::InfectionStatus = 0

num_infected

positionCoordinate x

pathMain.workspace3.Coordinate1Coordinate x

2Coordinate x

2Agent status

S::InfectionStatus = 0

num_infected

positionCoordinate x

pathMain.workspace3.Coordinate1Coordinate x

2Coordinate x

3Coordinate x

3Agent status

S::InfectionStatus = 0

num_infected

positionCoordinate x

pathMain.workspace3.Coordinate1Coordinate x

4Agent status

S::InfectionStatus = 0

num_infected

positionCoordinate x

pathMain.workspace3.Coordinate1Coordinate x

2Coordinate x

5Agent status

S::InfectionStatus = 0

num_infected

positionCoordinate x

pathMain.workspace3.Coordinate1Coordinate x

2Coordinate x

6Agent status

S::InfectionStatus = 0

num_infected

positionCoordinate x

pathMain.workspace3.Coordinate1Coordinate x

2Coordinate x

x
 
let # test step! function with toy example
    infection = CollisionInfectionRecovery(1., 0.002)
    agents = [Agent(S, 0, Coordinate(0,0)), Agent(S, 0, Coordinate(1,1)), Agent(S, 0, Coordinate(1,1)), Agent(S, 0, Coordinate(1,1)), Agent(S, 0, Coordinate(1,1)), Agent(S, 0, Coordinate(1,1))]
    for i in 1:length(agents)
        step!(agents, 5, infection)
        @debug "$(i): $(agents)"
    end
    agents
end

6.1 ms

Exercise 3.2

If we call step! N times, then every agent will have made one step, on average. Let's call this one sweep of the simulation.

👉 Create a before-and-after plot of $k_{s w e e p s} = 1000$ sweeps.

Initialize a new vector of agents (N=50, L=40, infection is given as pandemic below).
Plot the state using visualize, and save the plot as a variable plot_before.
Run k_sweeps sweeps.
Plot the state again, and store as plot_after.
Combine the two plots into a single figure using

plot(plot_before, plot_after)

20.3 μs

pandemic

CollisionInfectionRecovery p_infection

0.5

p_recovery

1.0e-5

6.3 ms

xxxxxxxxxx
 
@bind k_sweeps Slider(1:10000, default=1000)

85.1 ms

sweep! (generic function with 1 method)

xxxxxxxxxx
 
function sweep!(agent::Vector{Agent}, infection::AbstractInfection, L::Number)
    for i in 1:length(agent)
        step!(agent, L, infection)
    end
end

47.7 μs

xxxxxxxxxx
 
let
    N = 50  # number of agents
    L = 40  # plot size boundaries
    agents = initialize(N,L) 
    infection = pandemic
    plot_before = visualize(agents, L)
    for k in 1:k_sweeps
        # 1 sweep = N steps, where N = number of agents
        sweep!(agents, infection, L)
    end
    plot_after = visualize(agents, L)   
    plot(plot_before, plot_after)
 end

1.3 s

Exercise 3.3

Every time that you move the slider, a completely new simulation is created an run. This makes it hard to view the progress of a single simulation over time. So in this exercise, we we look at a single simulation, and plot the S, I and R curves.

👉 Plot the SIR curves of a single simulation, with the same parameters as in the previous exercise. Use k_sweep_max = 10000 as the total number of sweeps.

9.1 μs

k_sweep_max

xxxxxxxxxx
 
k_sweep_max = 10000

2.8 μs

simulation (generic function with 1 method)

x
 
function simulation(N::Integer, L::Integer, T::Integer, infection::AbstractInfection)
    agents = initialize(N, L)
    # track S I R after each sweep
    results = (S=[], I=[], R=[], R0=[0.])
    for k in 1:T
        sweep!(agents, infection, L)
        # tally S, I, R
        push!(results.S, count([a.status == S  for a in agents]))
        push!(results.I, count([a.status == I  for a in agents]))
        push!(results.R, count([a.status == R  for a in agents]))
        # calculate R0
        if k>1
            if results.I[k] > 0
                push!(results.R0, results.I[k] / results.I[k-1])
            else
                push!(results.R0, 0)
            end
        end
    end
    return results
end

144 μs

 
let
    N=40
    results = simulation(N, 30, k_sweep_max, pandemic)
    left = plot(1:k_sweep_max, results.S, ylim=N, label="S")
    plot!(left, results.I, ylim=N, label="I")
    plot!(left, results.R, ylim=N, label="R")
    right = plot(1:k_sweep_max, results.R0, label="R0")
    plot(left, right)
end

11.0 s

repeat_simulations (generic function with 1 method)

xxxxxxxxxx
 
function repeat_simulations(N::Integer, L::Integer, T::Integer, infection::AbstractInfection, num_simulations::Integer)
    map(1:num_simulations) do _
        simulation(N, L, T, infection)
    end
end

129 μs

plot_repeated_sims (generic function with 1 method)

xxxxxxxxxx
 
function plot_repeated_sims(simulation)
    n_time, n_sims = size(first(simulation).S)..., size(simulation)...
    left, right = plot(xlabel="time", ylabel="Number of agents Infected"), plot(xlabel="time", ylabel="R0")
    for sim in simulation
        plot!(left, 1:n_time, sim.I, alpha=0.5, label=nothing)
        plot!(right, 1:n_time, sim.R0, alpha=0.5, label=nothing)
    end
    i_mean = mean(reduce(hcat,[sim.I for sim in simulation]), dims=2)
    r0_mean = mean(reduce(hcat,[sim.R0 for sim in simulation]), dims=2)
    plot!(left, i_mean, lw=3, label=nothing)
    plot!(right, r0_mean, lw=3, label=nothing)
    plot(left,right, layout=(2,1))
end
​

159 μs

plot_repeated_sims2 (generic function with 1 method)

xxxxxxxxxx
 
function plot_repeated_sims2(simulation)
    n_time, n_sims = size(first(simulation).S)..., size(simulation)...
    p = plot(xlabel="time", ylabel="Number of agents")
    for sim in simulation
        plot!(p, 1:n_time, sim.I, alpha=0.5, label=nothing)
    end
    s_mean = mean(reduce(hcat,[sim.S for sim in simulation]), dims=2)
    i_mean = mean(reduce(hcat,[sim.I for sim in simulation]), dims=2)
    r_mean = mean(reduce(hcat,[sim.R for sim in simulation]), dims=2)
    plot!(p, s_mean, lw=3, label="S")
    plot!(p, i_mean, lw=3, label="I")
    plot!(p, r_mean, lw=3, label="R")
end
​

168 μs

Hint

After every sweep, count the values $S$ , $I$ and $R$ and push! them to 3 arrays.

150 ms

Exercise 3.4 (optional)

Let's make our plot come alive! There are two options to make our visualization dynamic:

👉1️⃣ Precompute one simulation run and save its intermediate states using deepcopy. You can then write an interactive visualization that shows both the state at time $t$ (using visualize) and the history of $S$ , $I$ and $R$ from time $0$ up to time $t$ . $t$ is controlled by a slider.

👉2️⃣ Use @gif from Plots.jl to turn a sequence of plots into an animation. Be careful to skip about 50 sweeps between each animation frame, otherwise the GIF becomes too large.

This an optional exercise, and our solution to 2️⃣ is given below.

9.5 μs

x
 
let
    N = 50
    L = 40
​
    agents = initialize(N, L)
    
    # initialize to empty arrays
    Ss, Is, Rs = Int[], Int[], Int[]
    
    Tmax = 200
    
    @gif for t in 1:Tmax
        for i in 1:50N
            step!(agents, L, pandemic)
        end
​
        #... track S, I, R in Ss Is and Rs
        sweep_result = StatsBase.countmap([agent.status for agent in agents])
        push!(Ss, get(sweep_result, S, 0))
        push!(Is, get(sweep_result, I, 0))
        push!(Rs, get(sweep_result, R, 0))
        
        left = visualize(agents, L)
    
        right = plot(xlim=(1,Tmax), ylim=(1,N), size=(600,300))
        plot!(right, 1:t, Ss, color=color(S), label="S")
        plot!(right, 1:t, Is, color=color(I), label="I")
        plot!(right, 1:t, Rs, color=color(R), label="R")
    
        plot(left, right)
    end
end

47.1 s

xxxxxxxxxx
 
let 
    N = 50
    for i in 1:50N
        println(i)
    end
end
​

77.9 ms

Hint

let
    N = 50
    L = 40

    x = initialize(N, L)
    
    # initialize to empty arrays
    Ss, Is, Rs = Int[], Int[], Int[]
    
    Tmax = 200
    
    @gif for t in 1:Tmax
        for i in 1:50N
            step!(x, L, pandemic)
        end

        #... track S, I, R in Ss Is and Rs
        
        left = visualize(x, L)
    
        right = plot(xlim=(1,Tmax), ylim=(1,N), size=(600,300))
        plot!(right, 1:t, Ss, color=color(S), label="S")
        plot!(right, 1:t, Is, color=color(I), label="I")
        plot!(right, 1:t, Rs, color=color(R), label="R")
    
        plot(left, right)
    end
end

10.0 μs

Exercise 3.5

👉 Using $L = 20$ and $N = 100$ , experiment with the infection and recovery probabilities until you find an epidemic outbreak. (Take the recovery probability quite small.) Modify the two infections below to match your observations.

8.1 μs

p_recovery0

5.0e-5

xxxxxxxxxx
 
p_recovery0 = 0.00005

3.0 μs

0.15

xxxxxxxxxx
 
@bind p_infect0 Slider(0.0001:0.0001:0.75, default=0.15, show_value=true)

345 ms

xxxxxxxxxx
 
let #not pandemic
N, L, T, n_sims = 100, 20, 500, 20
plot_repeated_sims(repeat_simulations(N, L, T, CollisionInfectionRecovery(p_infect0, p_recovery0), n_sims))
end

30.4 s

0.55

xxxxxxxxxx
 
@bind p_infect1 Slider(0.0001:0.0001:0.75, default=0.55, show_value=true)

102 μs

xxxxxxxxxx
 
let #pandemic
N, L, T, n_sims = 100, 20, 500, 20
plot_repeated_sims(repeat_simulations(N, L, T, CollisionInfectionRecovery(p_infect1, p_recovery0), n_sims))
end

23.5 s

causes_outbreak

CollisionInfectionRecovery p_infection

0.15

p_recovery

5.0e-5

xxxxxxxxxx
 
causes_outbreak = CollisionInfectionRecovery(0.15, p_recovery0)

9.3 μs

does_not_cause_outbreak

CollisionInfectionRecovery p_infection

0.55

p_recovery

5.0e-5

xxxxxxxxxx
 
does_not_cause_outbreak = CollisionInfectionRecovery(0.55, p_recovery0)

9.5 μs

Exercise 3.6

👉 With the parameters of Exercise 3.2, run 50 simulations. Plot $S$ , $I$ and $R$ as a function of time for each of them (with transparency!). This should look qualitatively similar to what you saw in the previous homework. You probably need different p_infection and p_recovery values from last week. Why?

8.3 μs

xxxxxxxxxx
 
let #pandemic
N, L, T, n_sims = 50, 40, k_sweeps, 50
plot_repeated_sims2(repeat_simulations(N, L, T, pandemic, n_sims))
end

68.9 s

need_different_parameters_because

The nature in which agents are interacting is fundamentally different, meaning the transmission of the diseases is fundamentally different. In the previous homework, agents interacted based on drawing from a uniform propbability distribution. In this homework, agents are interacting based on a random walk within a bounded space

xxxxxxxxxx
 
need_different_parameters_because = md"""
The nature in which agents are interacting is fundamentally different, meaning the transmission of the diseases is fundamentally different. In the previous homework, agents interacted based on drawing from a uniform propbability distribution. In this homework, agents are interacting based on a random walk within a bounded space
"""

16.5 μs

995 ns

Exercise 4: Effect of socialization

In this exercise we'll modify the simple mixing model. Instead of a constant mixing probability, i.e. a constant probability that any pair of people interact on a given day, we will have a variable probability associated with each agent, modelling the fact that some people are more or less social or contagious than others.

13.1 μs

Exercise 4.1

We create a new agent type SocialAgent with fields position, status, num_infected, and social_score. The attribute social_score represents an agent's probability of interacting with any other agent in the population.

11.3 μs

SocialAgent

x
 
# struct SocialAgent here...
begin
    mutable struct SocialAgent <: AbstractAgent
        status::InfectionStatus
        num_infected::Int64
        position::Coordinate
        social_score::Float64
        path::Array  # coordinate history
    end
​
    SocialAgent() = SocialAgent(S, 0, Coordinate(0,0), rand(LinRange(0.1,0.5,10)), [Coordinate(0,0)]) # default starting position is origin
    SocialAgent(s::InfectionStatus, n::Int64, p::Coordinate) = SocialAgent(s,n,p, rand(LinRange(0.1,0.5,10)), [p]) # starting point will always be the first coordinate in path history
end

12.6 ms

define the position and color methods for SocialAgent as we did for Agent. This will allow the visualize function to work. on both kinds of Agents

18.6 μs

👉 Create a function initialize_social that takes N and L, and creates N agents within a 2L x 2L box, with social_scores chosen from 10 equally-spaced between 0.1 and 0.5. (see LinRange)

10.4 μs

initialize_social (generic function with 1 method)

 
 function initialize_social(N::Number, L::Number)
    agents = [SocialAgent(S, 0, 
                          Coordinate(rand(-L:L),rand(-L:L))) for i in 1:N]
    rand(agents).status = I
    return agents
 end

84.3 μs

Now that we have 2 agent types

let's create an AbstractAgent type
Go back in the notebook and make the agent types a subtype of AbstractAgent.

24.2 μs

Exercise 4.2

Not all two agents who end up in the same grid point may actually interact in an infectious way – they may just be passing by and do not create enough exposure for communicating the disease.

👉 Write a new interact! method on SocialAgent which adds together the social_scores for two agents and uses that as the probability that they interact in a risky way. Only if they interact in a risky way, the infection is transmitted with the usual probability.

25.7 μs

interact! (generic function with 2 methods)

x
function interact!(agent::SocialAgent, source::SocialAgent, infection::CollisionInfectionRecovery)
    # agents must be in the same location to interact
    if agent.position == source.position
        # model risky interaction
        if rand() <= (agent.social_score + source.social_score)
            # for potential infect, one must be S and the other I 
            if agent.status==S && source.status==I && rand() <= infection.p_infection
                agent.status = I
                source.num_infected += 1
            end
        end
    elseif (agent.status == I)  && (rand() <= infection.p_recovery)
            agent.status = R 
    end
end 
​
​

98.2 μs

Main.workspace52.SocialAgent1SocialAgent status

I::InfectionStatus = 1

num_infected

positionCoordinate x

-1

social_score

0.45555555555555555

pathMain.workspace3.Coordinate1Coordinate x

2Coordinate x

3Coordinate x

-1

4Coordinate x

5Coordinate x

-1

6Coordinate x

-1

7Coordinate x

-2

8Coordinate x

-2

9Coordinate x

-2

10Coordinate x

-2

11Coordinate x

-1

-2

12Coordinate x

-2

-2

13Coordinate x

-2

-2

14Coordinate x

-1

-2

15Coordinate x

-2

16Coordinate x

-2

17Coordinate x

-2

18Coordinate x

-1

19Coordinate x

-1

-1

20Coordinate x

-2

-1

21Coordinate x

-2

-1

22Coordinate x

-2

23Coordinate x

-2

-1

24Coordinate x

-1

-1

25Coordinate x

-1

2SocialAgent status

I::InfectionStatus = 1

num_infected

positionCoordinate x

-1

-2

social_score

0.32222222222222224

pathMain.workspace3.Coordinate1Coordinate x

2Coordinate x

3Coordinate x

4Coordinate x

5Coordinate x

6Coordinate x

7Coordinate x

-1

8Coordinate x

-1

9Coordinate x

-2

10Coordinate x

-1

11Coordinate x

-1

12Coordinate x

13Coordinate x

14Coordinate x

15Coordinate x

-1

16Coordinate x

17Coordinate x

-1

18Coordinate x

-1

19Coordinate x

-1

20Coordinate x

-1

-1

21Coordinate x

-1

-2

22Coordinate x

-1

-2

23Coordinate x

-2

-2

24Coordinate x

-2

-2

25Coordinate x

-1

-2

26Coordinate x

-1

-1

27Coordinate x

-1

-2

x
 
let # test step! function with toy example
    infection = CollisionInfectionRecovery(1., 0.002)
    agents = [SocialAgent(I, 0, Coordinate(0,0)), SocialAgent(S, 0, Coordinate(1,1))]
    for i in 1:50
        step!(agents, 2, infection)
        @debug "$(i): $(agents)"
    end
    agents
end

799 ms

SocialAgent status

S::InfectionStatus = 0

num_infected

positionCoordinate x

social_score

0.4111111111111111

pathMain.workspace3.Coordinate1Coordinate x

2Coordinate x

3Coordinate x

4Coordinate x

5Coordinate x

6Coordinate x

7Coordinate x

8Coordinate x

9Coordinate x

10Coordinate x

11Coordinate x

xxxxxxxxxx
 
begin  # test rand_move! function on SocialAgent
a = SocialAgent()
rand_move!(a, 10, 5)
a
end

175 μs

Coordinate x

xxxxxxxxxx
 
position(a)

7.6 μs

Make sure step!, position, color, work on the type SocialAgent. If step! takes an untyped first argument, it should work for both Agent and SocialAgent types without any changes. We actually only need to specialize interact! on SocialAgent.

Exercise 4.3

👉 Plot the SIR curves of the resulting simulation.

N = 50; L = 40; number of steps = 200

In each step call step! 50N times.

23.5 μs

x
 
let
    N = 50
    L = 40
​
    global social_agents = initialize_social(N, L)
    Ss, Is, Rs = [], [], []
    
    Tmax = 200  # time period Tmax
    
    @gif for t in 1:Tmax
        # 1. Step! a lot
        for i in 1:100N  # for each time period step 100 sweeps
            step!(social_agents, L, pandemic)
        end
        # 2. Count S, I and R, push them to Ss Is Rs
        results = StatsBase.countmap([a.status for a in social_agents])
        push!(Ss, get(results, S, 0))
        push!(Is, get(results, I, 0))
        push!(Rs, get(results, R,0))
        
        # 3. call visualize on the agents,
        left = visualize(social_agents, L)
        # 4. place the SIR plot next to visualize.
        right = plot(xlim=(1,Tmax), ylim=(1,N))
        plot!(right, 1:t, Ss, color=color(S), label="S")
        plot!(right, 1:t, Is, color=color(I), label="I")
        plot!(right, 1:t, Rs, color=color(R), label="R")
        
        plot(left, right, size=(600,300)) # final plot
    end
end

68.6 s

Exercise 4.4

👉 Make a scatter plot showing each agent's social_score on one axis, and the num_infected from the simulation in the other axis. Run this simulation several times and comment on the results.

8.1 μs

plot_social_scores (generic function with 1 method)

x
 
function plot_social_scores(sagents::Vector{SocialAgent})
    social_scores, num_infected = Float64[], Int[]
    for a in sagents
        push!(social_scores, a.social_score)
        push!(num_infected, a.num_infected)
    end
    scatter(social_scores, num_infected, xlabel="social score", ylabel="number infected")
end

114 μs

x
 
plot_social_scores(social_agents)

126 ms

x
 
plot_social_scores(lockdown_agents)

13.6 ms

👉 Run a simulation for 100 steps, and then apply a "lockdown" where every agent's social score gets multiplied by 0.25, and then run a second simulation which runs on that same population from there. What do you notice? How does changing this factor form 0.25 to other numbers affect things?

7.5 μs

47.5 s

Incorporating a lockdown factor of 0.25 completely eliminates transmission. Loosening the lockdown factor to 0.75 results in a small uptick in transmission but it avoids a pandemic outcome.

13.6 μs

Exercise 5: (Optional) Effect of distancing

We can use a variant of the above model to investigate the effect of the mis-named "social distancing" (we want people to be socially close, but physically distant).

In this variant, we separate out the two effects "infection" and "movement": an infected agent chooses a neighbouring site, and if it finds a susceptible there then it infects it with probability $p_{I}$ . For simplicity we can ignore recovery.

Separately, an agent chooses a neighbouring site to move to, and moves there with probability $p_{M}$ if the site is vacant. (Otherwise it stays where it is.)

When $p_{M} = 0$ , the agents cannot move, and hence are completely quarantined in their original locations.

👉 How does the disease spread in this case?

40.8 μs

Enter cell code...xxxxxxxxxx
 
​

640 ns

👉 Run the dynamics repeatedly, and plot the sites which become infected.

6.6 μs

Enter cell code...xxxxxxxxxx
 
​

635 ns

👉 How does this change as you increase the density $ρ = N / (L^{2})$ of agents? Start with a small density.

This is basically the site percolation model.

When we increase $p_{M}$ , we allow some local motion via random walks.

17.8 μs

Enter cell code...xxxxxxxxxx
 
​

755 ns

👉 Investigate how this leaky quarantine affects the infection dynamics with different densities.

9.4 μs

Enter cell code...xxxxxxxxxx
 
​

634 ns

1.4 μs

8.8 μs

Function library

Just some helper functions used in the notebook.

18.0 μs

hint (generic function with 1 method)

50.8 μs

almost (generic function with 1 method)

32.1 μs

still_missing (generic function with 2 methods)

97.0 μs

keep_working (generic function with 2 methods)

57.2 μs

yays

Markdown.MD1

Fantastic!

Splendid!

Great!

Yay ❤

Great! 🎉

Well done!

Keep it up!

Good job!

Awesome!

You got the right answer!

Let's move on to the next section.

78.6 ms

correct (generic function with 2 methods)

50.7 μs

not_defined (generic function with 1 method)

50.9 μs

UndefVarError: starting not defined

top-level scope@Local: 1

xxxxxxxxxx
 
trajectory(starting, steps, 30)

---

4.7 ms